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Abstract — This paper investigates the use of analytical algo- 
rithms to quantify the uncertainty in the remaining useful life 
(RUL) estimate of components used in aerospace applications. 
The prediction of RUL is affected by several sources of uncer- 
tainty and it is important to systematically quantify their com- 
bined effect by computing the uncertainty in the RUL prediction 
in order to aid risk assessment, risk mitigation, and decision- 
making. While sampling-based algorithms have been conven- 
tionally used for quantifying the uncertainty in RUL, analytical 
algorithms are computationally cheaper and sometimes, are 
better suited for online decision-making. While exact analytical 
algorithms are available only for certain special cases (for e.g., 
linear models with Gaussian variables), effective approxima- 
tions can be made using the the first-order second moment 
method (FOSM), the first-order reliability method (FORM), and 
the inverse first-order reliability method (Inverse FORM). These 
methods can be used not only to calculate the entire probability 
distribution of RUL but also to obtain probability bounds on 
RUL. This paper explains these three methods in detail and 
illustrates them using the state-space model of a lithium-ion 
battery. 
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1. Introduction 

The need for an accurate and efficient health management 
system has become exceedingly important in safety-critical 
and mission-critical aerospace systems. The most impor- 
tant goal of health management is to constantly monitor 
the performance of these aerospace systems, identify faults 
(diagnosis), predict possible failures in the near future, and 
quantify the remaining useful life (prognosis) in order to aid 
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online decision-making. Sometimes, it may be challenging 
to perform health monitoring on the whole system due to 
its sheer complexity, and therefore, diagnosis and prognosis 
need to be performed on individual components that con- 
stitute the overall system. In this approach, mathematical 
models are developed for individual components, and then 
the component models are integrated to form the overall 
system. These component-level mathematical models can 
be constructed either using laws of physics (physics-based 
models [1]) or using data collected through component-level 
testing (data driven models [2]), and are used for both system- 
level diagnostics [3] and prognostics [4]. 

Uncertainty management is an important aspect of health 
monitoring, due to the presence of several unknown factors 
that affect the operations of the system of interest. Therefore, 
it is not only important to develop robust algorithms for 
diagnosis and prognosis, i.e., accurately perform diagnosis 
and prognosis in the presence of uncertainty, but also im- 
portant to quantify the amount of confidence in the results 
of diagnosis and prognosis. This can be accomplished by 
quantifying the uncertainty in fault diagnosis and prognosis 
(future performance prediction and remaining useful life). It 
is important to perform such uncertainty quantification (UQ) 
online so as to enable in-flight decision-making capabilities. 
Sankararaman and Mahadevan [5, 6] developed statistical 
(both frequentist and Bayesian) approaches to quantify the 
uncertainty in the three steps of diagnosis (detection, isola- 
tion, estimation) in an online health monitoring framework. 
There have also been a few papers [7—10] which discuss 
uncertainty propagation in prognosis; however, many of these 
papers are either suitable only for offline prognosis or they do 
not provide a comprehensive treatment of uncertainty. For 
example, the “Damage Prognosis Project” at Los Alamos 
National Laboratory [9] dealt with prognosis mainly in the 
context of offline testing and decision-making. Guan et 
al. [1 1 1 investigated Bayesian and maximum relative entropy 
methods for continuously updating the uncertainty in damage 
assessment, in the context of offline prognosis. Often, all 
the different sources of uncertainty - physical variability, 
data uncertainty, and model uncertainty - are not rigorously 
accounted for during prognosis; while most studies focus 
only on parameter variability and loading variability, the other 
sources of uncertainty are ignored. Sankararaman et al. [7] 
explained and addressed this issue in detail by identifying 
and accounting for the different types of uncertainty, the 
methodology was still developed in the context of offline 
testing. Therefore, there are several challenges relating to the 
topic of prognostics uncertainty quantification, and it is clear 
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that further research is necessary in this regard. The rest of the 
paper focuses on this aspect, i.e., quantifying the uncertainty 
in the reamining useful life (RUL) prediction. 

This paper pursues a model-based framework for prognosis 
and RUL calculation. Since state-space models are suitable 
for representing time-dependent behavior, the governing dif- 
ferential equations of the component of interest are trans- 
formed into equivalent state space representations. There are 
three major challenges before prognosis can be performed 
and the RUL can be calculated. 

1 . State Uncertainty: The first challenge is to estimate the 
current state of the system and the associated uncertainty. 
Typically, this problem is formulated as an estimation prob- 
lem or a state identification problem, and filtering techniques 
such as particle filtering [12] are used. Though the true state 
at any time instant is a deterministic quantity, i.e., it is not 
precisely known and in most cases, impossible to estimate 
with certainty due to (i) the presence of measurement errors 
and sensor noise (data uncertainty) in the data collected for 
estimation; (ii) uncertainty in the initial state; (iii) the noise 
in the state space model used for estimation; and (iv) model 
uncertainty. Filtering techniques such as particle filtering are 
based on the Bayesian philosophy of uncertainty representa- 
tion, according to which randomness can be perceived as lack 
of information/precision [13], and can be used to represent 
even an inherently deterministic quantity (which does not 
exhibit natural variabililty but is not known precisely) using a 
probability distribution. 

2. Future Loading Uncertainty: In practical engineering 
systems, it is almost impossible to accurately predict the 
future loading and environmental conditions. Following the 
Bayesian approach, future loading needs to be represented as 
an uncertain quantity. An important challenge is to quantify 
the amount of uncertainty in future loading. Though a 
few publications have addressed this issue [14, 15], further 
research is needed in this direction. However, the focus of 
the present paper is not on uncertainty characterization, and 
therefore it is assumed that future loading variability has been 
characterized (e.g., based on already existing loading profiles 
and data sets). After this uncertainty is characterized and 
quantified, it needs to be included in prognosis to quantify 
its effect on the uncertainty in RUL. 

3. Process Noise: Process noise is an important component 
of model uncertainty, and commonly represented as a random 
variable which needs to be included while future predictions 
are made and the RUL is estimated. Though the process noise 
term can be estimated using the data collected during the 
estimation stage, the future process noise may be significantly 
different from the estimated noise, and cannot be known in 
advance. However, it is still necessary to assume statistical 
distributions for the process noise in the different states, 
and include these distributions in the calculation of RUL 
uncertainty. 


The task of quantifying the uncertainty in the RUL prediction 
is primarily an uncertainty propagation problem [16], i.e., 
the above discussed different types of uncertainty need to be 
propagated through the state space model until the End of Life 
(EOL). Sampling-based methods (adaptive sampling, Latin 
hypercube sampling, Markov Chain Monte Carlo sampling, 
etc.) can be easily used to quantify the uncertainty in RUL, 
but they have two major disadvantages: 


1. Several thousands of samples and hence, several thou- 
sands of prognostic evaluations may be necessary to accu- 
rately quantify the entire probability distribution of RUL, 
which may not be computationally feasible for online health 
monitoring. Theoretically, an infinite number of samples is 
necessary to accurately calculate the probability distribution 
of RUL. An increasingly smaller set of samples result in 
increasingly worse representation of the distribution. 

2. The method of drawing random samples may not be 
preferred if repeatability is a desired criterion (for example, 
for online decision-making or also for certification of algo- 
rithms). In other words, the PDF of RUL depends on the 
exact set of considered samples; if the calculation is repeated 
with a different set of samples (which is likely the case if 
purely random samples are drawn), then a different (albeit 
only slightly, if the number of samples is large) probability 
distribution may be obtained. This difficulty can be overcome 
by predetermining the percentile values which can in turn be 
used to draw samples. In this case, it may be possible to 
preserve a particular statistic of the probability distribution, 
but it may be difficult to accurately obtain the entire probabil- 
ity distribution. For example, using the unscented transform 
sampling approach proposed by Daigle et al. [17], the mean 
and the variance of RUL can be efficiently calculated using 
a few samples. While it is important to calculate the first 
two moments of RUL, it is equally important to estimate the 
tail probabilities of the distribution of RUL, because failure 
is often caused due to events related to the tails of probability 
distributions. In other words, in a well-designed system, the 
probability of failure is very small (usually of the order of 
10 -3 or less), and in order to compute the value of RUL 
which corresponds to such low probabilities, it is important 
to accurately estimate the tail of the probability distribution 
of RUL. 


This paper investigates the use of analytical algorithms for 
calculating the uncertainty in RUL, as an alternative to 
sampling-based methods. Several analytical algorithms such 
as the first-order second moment method (FOSM), first-order 
reliability method (FORM), second-order reliability method 
(SORM), etc. have been used to calculate the reliability of 
structural systems, and they can also be used for uncertainty 
propagation. This paper investigates the use of the FOSM 
and FORM methods using state space models, and extends 
them to calculate the uncertainty in RUL prediction. These 
methods not only require very few prognostic evaluations in 
comparison with the sampling-based approaches but can also 
produce repeatable calculations, i.e., the exact same PDF on 
every repetition of the algorithm. While the former directly 
aids in online prognosis since fewer evaluations would lead to 
quicker calculations, it is worth noting that the latter feature is 
an important criterion for existing verification, validation, and 
certification protocols in the aerospace domain. Therefore, 
investigating such analytical algorithms allows us to move a 
step closer towards adopting prediction algorithms (which are 
inherently stochastic), by meeting the needs of the current 
certification process. 

The rest of the paper is organized as follows. Section 2 
explains the state space model formulation and describes how 
the estimation of uncertainty in RUL can be viewed as an 
uncertainty propagation problem. Sections 3 5 present the 
different analytical methods for quantifying the uncertainty 
in the RUL prediction. Finally, the methods are illustrated 
using the state space model of a lithium-ion battery model in 
Section 6, and Section 7 concludes the paper. 
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2. Prognosis and RUL Calculation 

The section formulates the prognosis problem, and explains 
the calculation of Remaining Useful Life [4], The various 
sources of uncertainty are explained, and finally, it is il- 
lustrated how the calculation of uncertainty in RUL can be 
treated as an uncertainty propagation problem. 

To begin with, consider the state space model which is used 
to continuously predict the state of the system, and therefore, 
aid in prognosis. Consider a generalized state space model 
as: 

x(t) = f (t, x(t), 0(t),u(t),v(t)) (1) 

where x(£) £ R"* is the state vector, 0(t) £ R n ® is the 
parameter vector, u(£) £ R"" is the input vector, v(f ) £ R ra ” 
is the process noise vector, f is the state equation. Typically, 
all of these quantities are uncertain in reality. 

The state vector at time £, i.e., x(t) (and the parameters 
0(t), if they are unknown) is (are) typically estimated using 
filtering approaches until the time t when data is available. 
Let y (£ ) £ R"“ , n (£) £ R n " , and h denote the output vector, 
measurement noise vector, and output equation respectively. 
Then, 

y(f) = h(£, x(£), 0(f), u(£), n(£)) (2) 

Note that the output equation or data is not used in the 
prognosis stage, since the focus is on predicting the future 
and the associated uncertainty. 


Prognostics and RUL prediction are concerned with the per- 
formance of the component that lies outside a given region 
of acceptable behavior. The desired performance is ex- 
pressed through a set of n c constraints, Ceol = 
where Ci : R"* x R rae x R”“ B maps a given point 
in the joint state-parameter space given the current inputs, 
(x(£), 0(t), u(£)), to the Boolean domain B = [0, 1], where 
Cj(x(f), 0(t), u(£)) = 1 if the constraint is satisfied, and 0 
otherwise. 


These individual constraints may be combined into a single 
threshold function Teol ■ R nx x R n ® x R n “ — > B, defined 


T E OL(x(t),0(t),u(t)) = I*’ 


0 G {ci(x(t),0(t), ujt))}"^ 
otherwise. 


( 3 ) 


Teol is equal to 1 when any of the constraints are violated. 
Then, the End of Life (EOL, denoted by E) at any time instant 
tp is then defined as the earliest time point at which this 
occurs: 

E(t P ) = inf{t G R : t > t P A T E oL( x (t), 0(t),u(t)) = 1} (4) 


The Remaining Useful Life (RUL, denoted by R ) at time 
instant tp is expressed as: 


R(tp) — E(t P ) — t P . (5) 


Thus, it is clear that RUL predicted at time tp, i.e., R(tp) 
depends on 

1. Present time (fp) 

2. Present state estimate (x(£p)); using the present state 
estimate and the state space equations in Eq. 1, the future 
states (x(fp), x(£p + 1), x(ip + 2), ..., x(ip + R(tp))) can 
be calculated. 


3. Future loading (u(£p), u(£p + 1), u(£p + 2), .... u(tp + 
R(tp))); these values are needed to calculate the future 
state values using the state space equations. In this paper, 
constant amplitude loading is considered, and the amplitude 
is assumed to follow a particular probability distribution. 
Future work will consider variable amplitude loading. 

4. Parameters (0(tp : tp + R(tp)))\ the values of parameters 
are also needed to calculate the future states. If measurement 
data are available, then the parameters can be estimated along 
with the states using filtering techniques. However, this paper 
addresses only the prognosis problem (future prediction), and 
it is assumed that measurement data is not available after time 
tp. Therefore, the parameters 9{tp : fp + R(tp )) need to be 
assumed in order to calculate future states and estimate the 
RUL prediction. In this paper, for the sake of illustration, 
the parameters are assumed to be constants (over time) and 
precisely known. 

5. Process noise (v(fp), v(£p + 1), v(£p + 2), ..., v(fp + 

R(t P ))y 

The focus of this paper is to quantify the uncertainty in RUL 
as a result of the uncertainties in these quantities. Note 
that the dependence of R(tp) on other uncertain quantities 
is implicit, i.e., the RUL depends on the process noise and 
loading, and in turn, how long the process noise and the fu- 
ture loading need to be considered during prognosis directly 
depends on the RUL itself. Therefore, it is challenging to 
employ the analytical methods for uncertainty propagation, 
since an explicit function, which establishes a one-to-one 
relationship from the various sources of uncertainty to the 
RUL, is necessary. The constant value of the input loading 
amplitude is denoted by up, and can be easily included as 
an input to this explicit function. However, the inclusion of 
process noise is not straightforward because it is necessary 
to include the process noise at all time instants as inputs; 
this increases the dimensionality of the problem several fold. 
Further, in some cases where the sensitivity of the uncertainty 
in RUL to the uncertainty in the process noise is negligible, 
it is meaningful to ignore the contribution of process noise 
uncertainty. Therefore, process noise is not considered in the 
present paper and the one-to-one function can be written as: 

R(t P ) = G(x(fp ), u E ) (6) 

Given a realization of x(fp) and u E , the above equation 
is continuously evaluated using the state space model in 
Eq. 1, and R(tp) is computed when end of life is reached. 
In a sampling-based procedure, this is repeated for several 
thousands of samples of x(fp) and u B , the corresponding 
realization of R(tp) are calculated to obtain a histogram, 
which can be used to construct the probability density func- 
tion (PDF) of R(tp). As stated earlier, such sampling-based 
approaches may be computationally expensive, and this paper 
investigates computationally efficient, analytical alternatives 
to construct the PDF. In general, there may be other sources 
of uncertainty (e.g., parameter uncertainty), and hence, the 
above function G may have several inputs. The following 
sections discuss the analytical methods for a generic number 
of inputs, as: 

R = G(X) (7) 

where X is a vector X = {A' 1 ,X 2 , ...Xj, ...X n }, where n 
is the length of the vector X, and therefore the number 
of uncertain “inputs” to G. For example, if there are two 
states (their values at time t p are uncertain) and one loading 
amplitude (that is uncertain), then n = 3. Each X, has its 
own probability distribution and it is desired to calculate the 
probability distribution of R. If these probability distributions 
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are Gaussian and if G is a linear function, then uncertainty 
propagation through G is simple and closed-form expressions 
are available to estimate the statistics of R. In general, 
the state-space equation, i.e., Eq. 1, and therefore, damage 
progression equations are non-linear; further, the threshold 
function in Eq. 3 may also be non-linear. Therefore, G is 
non-linear in several practical cases, and thus the uncertainty 
quantification methodology must be applicable to the general 
case, i.e, non-linearity in G and non-normality in X. 


2. FORM for Non-Normal Variables: The above FORM 
method for normal variables is extended to non-normal vari- 
ables in the latter part of Section 4. 

3. Inverse FORM: While Section 4 deals with calculating 
the CDF value for a given realization of R, Section 5 con- 
siders the inverse problem, i.e. calculating the realization 
of R that corresponds to a given CDF value. This method, 
popularly known as the Inverse CDF approach, can be used 
to calculate probability bounds on the remaining useful life 
prediction, which is useful for decision-making. 


3. First-Order Second Moment Method 


The first-order second moment approach, as the name sug- 
gests, is a simple approximation of R using First-order Taylor 
series expansion. The first two moments, i.e. mean (fix) 
and variance (ax) of X are used to approximate the first 
two moments of R. Consider the first-order Taylor series 
expansion of Y = G(X) around fix, as: 

R = G(fix) + (8) 

Note that R is a linear function of X with the partial deriva- 
tives as coefficients, and therefore, it is straightforward to 
approximate the mean and variance of R, as: 

HR = G(fix) ( 9 ) 


°R 


i=nj=n , „ 
*= 1 3 = 1 1 * * * V 


MX 


d R 

dX-i 


Cov(X l ,X j ) (10) 


MX 


When the inputs to G are uncorrelated, then the expression 
for variance in Eq. 10 simplifies to: 
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<Tx t ( 11 ) 

MX 


Note that the first-order second moment approach can only 
compute the mean and variance of the remaining useful life 
(R) prediction. It is not directly useful in estimating the 
type of probability distribution of R, and therefore cannot be 
used to calculate tail probabilities. However, if each of the 
inputs Xi were to be normally distributed (and independent 
of each other), then it can be easily proved that their linear 
combination is also a normally distributed variable [18]. This 
observation is used to develop a general methodology to 
calculate the entire CDF of the Remaining Useful Life (R) 
prediction in the following stages: 

1. FORM for Normal Variables: First, consider the case 

when the variables are normal. Even in this case, the above 

method is suitable only when G is linear. When G is non- 

linear, it is not appropriate to always calculate the gradient 

at the mean fix ■ This leads to obvious question — “where 
to linearize”? The “location of linearization” is chosen 
analytically, and then used to calculate the CDF of RUL at 
a particular R = r, i.e., Fr(R = r) = P(R < r). Note 
that an upper case letter refers to the name of a random 
variable whereas the corresponding lower case letter refers 
to its realization. It will be illustrated that the location (point) 
of linearization varies with the choice of r. By repeating the 
entire process for different choices of r, the entire CDF can be 
calculated. This method, known as the first-order reliability 
method (FORM), is discussed in the first part of Section 4. 


4. First-Order Reliability Method 

The first-order reliability method (FORM) was originally 
developed by structural engineers in order to estimate the 
reliability of structural systems. The contribution of this 
paper is to extend this method to state space models and 
health monitoring with the goal of computing the probability 
distribution of Remaining Useful Life. Though the estimation 
of reliability and the calculation of CDF are statistically 
equivalent, the FORM method is described from a purely 
uncertainty propagation perspective, without reference to re- 
liability calculation. 

Consider Eq. 7, which expresses the RUL as a function of the 
various sources of uncertainty. The first subsection explains 
the FORM method when X consists of normally distributed 
variables and the next subsection extends the method to the 
non-normal case. 

Normal Variables 

The goal of the FORM approach is to calculate the CDF 
value, i.e. Fn(r) = P(R < r), given r which is a realization 
of the random variable R. In this subsection, assume that the 
i th input Xi is normally distributed as N(fix,, ax,)- 

FORM achieves the aforementioned goal by approximating 
the non-linear equation R = G(X) using a linear equation in 
order to easily compute the CDF of R. The linear equation 
is constructed using a Taylor series approximation around the 
so-called “point of linearization”. The difference in FORM 
(with respect to FOSM) is that the point of linearization varies 
from one choice of r to another. The identification of the 
point of linearization is the most important component of the 
FORM algorithm. 

In order to calculate F R (r) = P(R < r), consider the 
contours of the function R — r = G{X) — r; in particular 
consider the curve described by G(x) — r = 0, where x 
is a realization of X. This curve differentiates the multidi- 
mensional space into two regions: one region where R < r, 
and another region where R > r, as shown in Fig. 1. (This 
curve of demarcation is popularly called as the limit state in 
reliability analysis. The equation corresponding to the limit 
state divides the space into two zones - zone of failure and 
zone of safety, respectively. This terminology is not used 
in this paper because FORM is now used for uncertainty 
propagation rather than reliability calculation.) 

Any point lying on the curve of demarcation would satisfy 
the equation R = G (X). Since (1) this curve serves as the 
demarcation between the two zones given by R > r and 
R < r; and (2) it is of interest to calculate the probability 
P(R < r), it is intuitive that it is important to identify a linear 
function which closely resembles the contour G(x) r = 0. 
Hence, the point of linearization must he on this curve of 
demarcation; in other words, the point of linearization must 
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Figure 1: Estimating MPP in FORM 


satisfy the equation G{x) — r = 0. This is clearly different 
from the FOSM approach, where the mean px was chosen 
as the point of linearization; for an arbitrary choice of r, it is 
obvious that the mean p x will not satisfy this equation. 


Eq. 14 can rewritten as: 

i—n 

= (15) 

4=1 

where 


<Ji 

If the above computation were performed for every realiza- 
tion Xi of the random variable X,, then the corresponding 
UiS would be realizations of the standard normal variable 
Uj , i.e., U, ~ iV(0, 1). Therefore, Eq. 16 is referred to as 
the standard normal transformation. In the space of standard 
normal variables, maximizing the likelihood of occurrence is 
equivalent to minimizing Eq. 15, which implies that the point 
of linearization is that point on the curve of demarcation, 
whose distance (measured in the standard normal space) from 
the origin is minimum. Since the point of linearization has the 
maximum likelihood of occurrence, it is popularly known as 
the Maximum Probable Point (MPP), as indicated in Fig. 1 . 
The Maximum Probable Point, therefore, represents the min- 
imum distance (measured from the origin, in the standard 
normal space), and this minimum distance is of significance, 
because it can be proved that: 


Therefore, the point of linearization should be located on the 
curve of demarcation. However, there are infinite points that 
satisfy this criterion, and it is important to select the appro- 
priate one. Each of these infinite points has a likelihood of 
occurrence, and intuitively, the point of maximum likelihood 
is chosen as the point of linearization. This likelihood can 
be calculated using the probability density function of the 
underlying random variables. For a single normal random 
variable X with mean p and standard deviation a, the PDF is 
given by: 


fx(x\p,a) 


T\/2tt 


exp 


{x- p) 2 ' 
2o 2 


(12) 


For example, when p = 10 and a = 1, x = 10 is 1.65 times 
more likely to occur than x = 9. The maximum value of the 
likelihood function occurs at x = p\ therefore, the farther x is 
away from the mean p, the lower the likelihood of occurrence 
of x. As explained earlier, the mean cannot be chosen as the 
point of linearization since G{p) — r f 0. Therefore, if there 
is a single input variable, the point of linearization is chosen 
in such a way that it satisfies the equation G(x) — r = 0 and 
the value of | |tc — p\ \ is minimum. 


However, in a general uncertainty propagation problem, the 
input to G is a vector, i.e., X = {Xi,X 2 , ...Xi, ...X n }, and 
each Xi has its own mean p Xi and standard deviation a Xi ■ 
The objective is to identify the point of maximum likelihood, 
which can be calculated by maximizing the joint probability 
density function of all the input random variables. If the 
variables are independent, then the joint density function of 
X is expressed as: 


fx(x) 


i=n 

n 


i 


-=exp 

CTiV 27T 


(Xi 


) 2 1 

2°f J- 


(13) 


It can be verified (by taking the logarithm) that the maximizer 
of the above function simultaneously minimizes 


i=n „ 


(14) 


P(R <r) = $(— /3) (17) 

where <f>(.) represents the standard normal cumulative dis- 
tribution function. In order to prove Eq. 17, consider the 
first order Taylor’s series expansion of Z = G(X) — r by 
linearizing around the MPP. Then, Z is approximated to be a 
normal random variable; its CDF value measured at Z = 0 is 
exactly equal to P(R < r), and is calculated as: 

P(Z < 0) = P(R < r) = <&(— — ) (18) 


Using the Taylor series expansion in Eq. 8, it can be proved 
that [18]: 




Pz_ 

crz 


(19) 


Therefore, the problem of calculating the CDF reduces to 
identifying the MPP on the curve of demarcation. This can 
be posed as a constrained optimization problem, as follows: 

Minimize u T u 

U 

s.t. G(x ) = r 

U = {zti, U2-..Ui, (20) 

Ui = (Xi - Pi)/Ui 

(i = 1 to n) 

The above optimization problem can be solved using the 
Rackwitz-Fiessler [19] algorithm, an iterative procedure, as 
follows: 


1. Initialize counter j = 0 and start with an ini- 

tial guess for the Most probable point (MPP), i.e., 
xi = , a? 2 , ... x%, a column vector. 

2. Transform into standard normal space and calculate 
u° = {«,{, u J 2 , uj, using Eq. 16, a column vector. 

3. Compute the gradient vector in the standard normal space, 

i.e., a = {cti, a 2 , another column vector where 


dg dg da dg 

dui dxi dui dxi 


(21) 
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4. In the iterative procedure, the next point u 3+1 is calculated 
using a Newton-Raphson type equation, as: 

u j+1 = [a T u j - G(x j )] — — — (22) 

Ha'll | |o:| | 

5. Transform back into original space, i.e., compute x J + 1 , 
and continue starting from Step 3 until the iterative procedure 
converges. Using tolerance limits and 82 , convergence 
can be checked verified if the following two criteria are 
satisfied: (i) the point lies on the curve of demarcation, 

i.e., |G(a; J ) — r| < ch; and (2) the solution does not change 
between two iterations, i.e., \x ^ +1 — xJ \ < 82 - 

The above described iterative procedure usually converges 
within 4 or 5 iterations. In each iteration, the transformation 
to the standard normal space is straightforward only when the 
variables are originally normal. Therefore, the method needs 
to be modified to include non-normal variables, as explained 
in the next subsection. 

Extension to Non-normal Variables 

Now consider the case where the inputs X, (i = 1 to n) 
have arbitrary probability distributions given by their CDFs as 
Fxi = 1 to n). Now that Xi is not normally distributed, 
Eq. 16 cannot be used for standard normal transformation. 
Therefore, it is necessary to calculate Ui from a given x t 
meaningfully, so that m, represents a realization of the stan- 
dard normal variable. One simple transformation is based on 
probability integral transform concept, as: 

u i = $- 1 (F Xi (X i =x i )) (23) 

where <b _1 (.) refers to the inverse of the standard normal 
distribution function [18]. 

In addition to the above procedure, there are also other 
transformation techniques. For example, a two-parameter 
transformation procedure estimates the mean fii and standard 
deviation a, of the normal distribution by equating the CDF 
and PDF values of the distribution of X and the normal distri- 
bution. Then, Eq. 16 can be used to calculate Ui from x^. Note 
that the mean //., and standard deviation < 7 i are dependent on 
the value of Xi. Similarly, Chen and Lind [20] proposed a 
three-parameter transformation procedure by introducing a 
third parameter, a scale factor which is estimated by matching 
the slope of the probability density function in addition to the 
PDF and CDF values. Further, when the inputs are corre- 
lated or statistically dependent, it is necessary to transform 
them to uncorrelated standard normal space. Haidar and 
Mahadevan [18] describe methods for such transformation. It 
must be noted that any transformation must be accompanied 
by suitably computing the derivatives in the standard normal 
space, and Eq. 21 must be appropriately replaced. 

Note that above FORM procedure calculates the CDF value 
at a particular value of RUL. It answers the question: What is 
the probability that the RUL is smaller than a given number? 
This question is answered using a search procedure which 
computes the CDf value corresponding to a given RUL value, 
and this search is accomplished through an iterative proce- 
dure, like in an optimization algorithm. In order to obtain the 
entire CDF, the whole procedure is repeated with multiple 
values of RUL. Sometimes, it may not be possible to identify 
values of RUL in order to calculate the entire CDF because 
the spread of the distribution may not be known in advance. 
So, the next section discusses the Inverse FORM procedure. 


which answers the question: What is the value of RUL which 
corresponds to a given probability level? In other words, 
what is the o-percentile (e.g., 5%, 95%, etc.) value of RUL? 
By repeating this procedure for one lower percentile and one 
upper percentile value, the probability bounds on RUL can be 
calculated. 


5. Inverse First-Order Reliability 
Method 

Given j3 or A = <&(— 0), the inverse FORM approach can be 
used to calculate r such that Fn{r) = P(R < r) = A. The 
theory behind inverse FORM is exactly the same as FORM, 
and the algorithm discussed in Section 4 is modified so that 
the CDF value can be specified and r can be calculated. The 
various steps involved in the iterative procedure for Inverse 
FORM are outlined below: 

1. Initialize counter j = 0 and start with an ini- 

tial guess for the Most probable point (MPP), i.e., 

3 / nr**) rpj "1 

0, — \«^i ? ^ 2 > J * 

2. Transform into standard normal space and calculate 

= {u{,u 3 2 ,... uj, 

3. Compute the gradient vector in the standard normal space, 

1. e., a = {ai, « 2 j as explained in Section 4. 

4. In the iterative procedure, the next point u :l+ i is calculated 
as: 

u 3+1 = -pr/3 (24) 


5. Transform back into original space, i.e., compute a: J + 1 , 
and continue starting from Step 3 until the iterative procedure 
converges. Using tolerance limits 8 \ and 82 , convergence 
can be checked verified if the following two criteria are 
satisfied: (i) the point lies on the curve of demarcation, 
i.e., | G{x 3 ) — r\ < <5i; and (2) the solution does not change 
between two iterations, i.e., \xJ +1 — xJ j < 82 - 

Similar to the Rackwitz-Fiessler algorithm, the above itera- 
tive procedure usually convergences within 4 or 5 iterations, 
and therefore is suitable for quick calculations. 

As stated earlier, the Inverse FORM procedure is useful 
to calculate probability bounds. For example, by repeat- 
ing the above algorithm for A = 0.05 and A = 0.95, 
it is possible to estimate the 90% probability bounds on 
the remaining useful life. In fact, the entire CDF can be 
constructed by repeating the analysis for several values of A 
(0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95). 

The advantage of Inverse FORM is that the number of prog- 
nostic evaluations is extremely small compared to sampling- 
based approaches. For the sake of illustration, consider that 
the above iterative algorithm converges in 5 iterations. Each 
iteration requires the computation of the gradient vector; if 
there are k inputs, then k + 1 computations are necessary 
for each evaluation - one for each derivative and one for 
the evaluation of G. Hence, the total number of prognostic 
evaluations is equal to 5 [k + 1). If the analysis is repeated for 
2 values of A (in order to compute the probability bounds on 
RUL), then it takes 1 ()(/,: + 1) prognostic evaluations, which 
is more computationally efficient than sampling-based tech- 
niques, and, therefore suitable for online decision-making. 
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6. Numerical Example 

The proposed methods are illustrated using a state space 
model of a lithium-ion battery, which is an important compo- 
nent in aerospace applications and which is also being studied 
at NASA Ames Research Center as part of the rover test- 
bed [21]. Though the method is illustrated using a battery 
model, it is general and can be applicable to state space 
models in several engineering domains. 

Description of the Model 

This battery model, earlier used by Daigle et al. [17] for 
prognosis, is similar to the models presented by Chen and 
Rincon-Mora [22]. The model is based on an electrical circuit 
equivalent as shown in Fig. 2, where the large capacitance 
C b holds the charge qt, of the battery. The Rcp-Ccp pair 
captures the major nonlinear voltage drop due to concentra- 
tion polarization, R s captures the so-called 1-R drop, and 
R p models the parasitic resistance that accounts for self- 
discharge. This relatively simple battery model is sufficient 
to capture the major dynamics of the battery, but ignores 
temperature effects and other minor battery processes. 


Rcp 



The state-of-charge, SOC, is computed as 

SOC = 1 - 9m ™ ~~ qb (25) 

{-'max 

where qb is the current charge in the battery (related to 
Cb), q m ax is the maximum possible charge, and C max is 
the maximum possible battery capacity (i.e., nominally, its 
rated capacity). The concentration polarization resistance is a 
nonlinear function of SOC: 

Rcp = Rcpo + Rcpi exp(i?cP2(l — SOC)) (26) 

where Rcpo , Rcpi , and Rcpi are empirical parameters. 
The resistance, and, hence, the voltage drop, increases ex- 
ponentially as SOC decreases [14]. 

Voltage drops across the individual circuit elements are given 
by 

V b = q b /Cb (27) 

Vcp = qcp/ Ccp (28) 

V p = V b - Vcp (29) 

where qcp is the charge associated with the capacitance 
Ccp ■ The terminal voltage of the battery is 

V = V b -V C p- R s i (30) 


Parameter 

Value 

Ub 

9844 

Rs 

0.143014 

R P 

500 

Ccp 

70.3767 

Rcpo 

0.019829 

Rcpi 

3.68606 x 10“ 14 

Rcp2 

31.9213 

Qmax 

41400 

Cmax 

6900 


Table 1: Battery Model Parameters 


where i is the battery current at the terminals. Currents 
associated with the individual circuit elements are given by 


i p = Vp/Rp (31) 

i b = ip + i (32) 

icp = ib — Vcp /Rcp (33) 

The charges are then governed by 

qb = —ib (34) 

qcp = icp (35) 


It is of interest to predict the end-of-discharge as defined by 
a voltage threshold Veod- So, Ceol consists of only one 
constraint: 

ci : V > Veod (36) 


The parameters of the battery model are assumed to be 
deterministic and are shown in Table 1. All voltages are 
measured in Volts, resistances are measured in Ohms, charges 
are in the units of Coulombs, and capacitances are measured 
in Coulombs per Volt. 

The following sections deal with uncertainty quantification 
and verification of the Inverse FORM approach. 

Uncertainty Quantification in RUL 

In this example, the process noise of the state space model 
is assumed to be small, and therefore not included in the 
uncertainty quantification. The two types of uncertainty 
considered are: 

1. Loading Uncertainty: In this example, a constant am- 
plitude loading is considered for the purpose of illustra- 
tion. The amplitude is considered to be normally distributed 
(_/V(1.375, 1/6)), and this distribution is truncated at a speci- 
fied lower bound (0.75) and an upper bound (2.00). 

2. State Uncertainty: Typically, the state estimation, which 
is an inverse problem, is addressed using a filtering technique 
which can continuously estimate the uncertainty in the state 
when measurements are continuously available as a function 
of time. In this paper, the state estimation is not explicitly 
carried out. The state values are assumed to be available, 
and the uncertainty in the states are predetermined based 
on the authors’ past experiences with the use of filtering 
techniques for the above described problem. There are two 
state variables (charge in the battery and charge associated 
with the capacitance) in this example and at any time instant, 
they are assumed to be normally distributed with a specified 
mean; for example, the mean of the initial states are set as 
[4.14 x 10 4 , 0]. For the purpose of illustration, three different 
values of CoV (Coefficient of variation, defined as the ratio 
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CoV=lA CoV=2% CoV=3% 

T=0 

ri = 4065 n = 3781 r i = 3401 

r 2 = 5045 r 2 = 5045 r 2 = 5045 

r 3 = 6459 r 3 = 6742 r 3 = 7120 

T=1000 

ri = 3216 n = 2908 ri = 2515 

r 2 = 4045 r 2 = 4045 r 2 = 4045 

r 3 = 5226 r 3 = 5532 r 3 = 5924 

T=2000 

r i = 2359 rj = 2023 ri = 1616 

r 2 = 3045 r 2 = 3045 r 2 = 3045 

r 3 = 4001 r 3 = 4335 r 3 = 4741 

T=3000 

ri = 1488 rj = 1118 ri = 699 

r 2 = 2045 r 2 = 2045 r 2 = 2045 

r 3 = 2789 r 3 = 3157 r 3 = 3575 

T=4000 

ri = 589 ri = 186 r\ = 43 

r 2 = 1045 r 2 = 1045 r 2 = 1045 

r 3 = 1605 r 3 = 2008 r 3 = 2432 

T=5000 

n r - 6 ri = 1 n — 0 

r 2 =45 r 2 = 45 r 2 = 45 

r 3 = 478 r 3 = 897 r 3 = 1316 


Table 2: Results of Uncertainty Quantification in RUL 



Figure 3: 90% Bounds and Median of RUL (CoV = 1%) 


between standard deviation and mean) are considered — 

0. 01, 0.02, 0.03 — and the analysis is repeated for each CoV 
value. 

In order to apply the proposed algorithm to this battery exam- 
ple, the first step is to construct the “G” model as seen earlier 
in Eq. 6. One realization of each source of uncertainty (two 
variables corresponding to the two states in this example, and 
one variable corresponding to the loading) becomes input 
to this model, and the state space equations above are used 
to calculate the remaining useful life prediction (right hand 
side of Eq. 6). Then, the Inverse FORM method can be 
directly used to compute the RUL value corresponding to 
multiple CDF values; Fr{t \ ) = 0.05, Fr(t 2 ) = 0.5, and 
Fr(t 3 ) = 0.95 are used in this paper, and the calculations are 
repeated for three considered CoV values for the uncertainty 
in the state variables. While n and r : > correspond to the 90% 
probability bounds of RUL, r 2 corresponds the median of 
RUL. The bounds and the mean are continuously calculated 
until T = 5000 (in seconds) when failure seems to be immi- 
nent, and the Inverse FORM calculation is performed every 
100 s. The results are tabulated in Table 2 and graphically 
shown in Fig. 3-5. 

It is seen that from the results that the uncertainty in the RUL 
is high initially, and then gradually decreases until failure is 
imminent. Initially, the uncertainty in RUL is high because it 
is necessary to predict at a further time instant, in comparison 
with making prognostic predictions at a later stage. In fact, 
any good prognostic algorithm should depict this behavior, 

1. e., the prediction of RUL at a later time instant has lower 



Figure 4: 90% Bounds and Median of RUL (CoV = 2%) 



Figure 5: 90% Bounds and Median of RUL: (CoV = 3%) 


uncertainty than the prediction at an earlier time instant. 
Further, the larger the coefficient of variation, the larger the 
uncertainty in RUL; this behavior is observed at every time 
instant and is consistent with intuition because a larger un- 
certainty in the state estimate implies that the corresponding 
RUL prediction uncertainty will also be high. When failure 
is imminent, the probability distribution of RUL is highly 
positively skewed; in other words, the probability density 
function has a very high likelihood for almost immediate 
failure, but has a longer right-tail. The maximum likelihood 
point corresponds to the time of imminent failure, while the 
variance is reflective of the various sources of uncertainty. 

Verification using Monte Carlo Sampling 

In order to verify the above performed uncertainty quantifica- 
tion, Monte Carlo sampling (MCS) is performed using 1000 
samples. It is computationally infeasible to perform MCS at 
every time instant considered. Therefore, MCS is performed 
for six time instants, starting from T = 0 until T = 5000 in 
steps of 1000 seconds. At each time instant, the CDF of RUL 
was computed by repeating Inverse FORM for 13 different A 
values (A = 0.01, 0.05, 0.1, 0.2, ...0.8, 0.9, 0.95, 0.99). The 
comparison of Inverse FORM and MCS at T = 0 is shown 
in Fig. 6; in this illustration, CoV of state estimates is chosen 
to be 3%. Note that the uncertainty bounds due to the use of 
limited number of samples for MCS are also shown. 

It is seen that the probability distribution resulting from 
Inverse FORM lies within the Monte Carlo bounds, as seen 
in Fig. 6, thus verifying the uncertainty quantification proce- 
dure. In fact, the maximum difference between the Inverse 
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RUL 


Figure 6: Inverse FORM vs. MCS (T = 0, CoV = 3%) 



Figure 7: Inverse FORM vs. MCS (T = 4000, CoV = 1%) 


FORM solution and the Monte carlo solution was found to 
be less than 0.5% (it must be noted that the Monte carlo 
solution is not exact due to the use of limited samples and 
hence, the 90% Monte Carlo bounds have also been provided 
in Fig. 6). Further, the Inverse FORM procedure needed 
only 13 x 4 x 3 = 156 prognostic evaluations, since 13 A 
values and 3 variables (to represent uncertain quantities) were 
used, and 4 iterations were needed for convergence whenever 
the iterative Inverse FORM algorithm was used. In contrast, 
MCS required 1000 prognostic evaluations. Further, a similar 
comparison was performed at different time steps, and by 
considering other CoV values. For instance, the comparison 
of the results from MCS and Inverse FORM at T — 4000 s 
for CoV = 1% is shown in Fig. 7. 

The agreement of Inverse FORM with MCS is evident in 
Fig. 7. In fact, it was observed that the Inverse FORM method 
performs well in comparison with Monte Carlo sampling, 
both in terms of accuracy (the maximum difference was found 
to be less than 1%) and computational cost (approximately 
10% of the cost of MCS). 


7. Conclusion 

Conventionally, sampling-based algorithms have been used 
for quantifying the uncertainty in prognostic calculations, and 
may require several thousands of samples in order to quantify 
the entire probability distribution of remaining useful life 
(RUL) with reasonable accuracy. Further, if such an algo- 
rithm were to be repeated, a probability distribution different 


(albeit only slightly) from the original distribution may be ob- 
tained, due to the use of limited samples. These two reasons 
hinder the use of such algorithms in online prognostics and 
decision-making. Therefore, this paper investigated the use 
of analytical algorithms in order to quantify the uncertainty 
in the RUL prediction. 

The First-Order Second Moment method (FOSM), the First- 
Order Reliability Method (FORM), and Inverse First-Order 
Reliability Method (Inverse FORM) are optimization-based 
analytical approaches, which were originally developed by 
structural engineers in order to quantify the reliability of 
structural systems. This paper extended these methods to 
quantify the entire probability distribution of RUL using 
state-space models that capture the evolution of the system 
continuously as a function of time. While the FOSM method 
can be used to obtain the first two moments of the RUL 
prediction, the FORM and Inverse FORM approaches can 
be used to calculate the entire CDF of the RUL. These 
approaches are not only computationally cheaper, but can 
produce repeatable results, and therefore, are suitable for 
online prognosis and decision-making. 

Future work needs to address several issues. First, the 
process noise of the state space models need to be rigorously 
accounted for in the uncertainty quantification procedure. 
This is particularly important when the contribution of pro- 
cess noise uncertainty to the overall uncertainty in RUL is 
significantly large. Second, practical systems are commonly 
subjected to different types of variable amplitude loading 
profiles such as block loading, Markov processes, general 
random processes, etc., and therefore, the proposed methods 
for uncertainty quantification need to be extended to consider 
variable amplitude loading. The assumption of constant am- 
plitude loading implies that loading uncertainty is described 
using a single random variable, whereas variable amplitude 
loading profiles need to be described using multiple random 
variables, which not only increases the dimensionality of 
the problem, but also affects the uncertainty in the RUL 
prediction. Future work must investigate the impact of in- 
cluding variable amplitude loading on the uncertainty bounds 
of the RUL prediction. Third, sensitivity analysis needs 
to be performed so that the contributions of the different 
sources of uncertainty to the overall uncertainty in RUL can 
be quantified. Fourth, this paper did not consider the effect of 
model form uncertainty on prognosis; future research needs 
to quantify model form uncertainty and develop a method to 
rigorously account for model uncertainty in prognosis and 
RUL calculations. Finally, it is also necessary to quantify the 
robustness of the proposed approach, by estimating the sensi- 
tivity of the RUL bounds to relaxing the various assumptions 
in the present paper, and thereby investigate the applicability 
of the methodology to practical engineering systems. 
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